% discrete legendre polynomials
clear; clc;
N = 15;      % max degree
k = 0:N;    % interval

P = zeros(N,N+1);
P(1,1:N+1) = 1; 
P(2,1:N+1) = 1 - (2*k)/N;
for i = 1:N-1
    P(i+2,1:N+1) = ( (2*i+1)*(N-2*k).*P(i+1,:) - i*(N+i+1)*P(i,:) ) ./ ((i+1)*(N-i)) ;
end

for i = 0:N % fading factorial:  FF(N,i) = N*(N-1)*...*(N-i+1) = N*(N-1)*... / (N-i)*(N-i-1)*... = F(N)/F(N-i);
    a = (2*i+1);
    b = factorial(N+i+1)/factorial(N);
    c = factorial(N)/factorial(N-i);
    w(i+1)=a*c/b;
end

% a function f(k) can be expanded:
%   f(k) = sum f(i)*p(i) = f'p
% where
%   f(i) = w(i) sum f(k)*p(i)
x = linspace(0,1,N+1);
f = x.*(-1).^(mod(k,2));
f = cos(x).*x;
% now we want to expand f into fi
ff(1) = w(1) * f * P(1,:)';
ff(2) = w(2) * f * P(2,:)';
ff(3) = w(3) * f * P(3,:)';
ff(4) = w(4) * f * P(4,:)';

% now approximate f with order 1
f1 = ff(1)*P(1,:);
plot(k,f,'k',k,f1,'r')
% and tada, it's the mean.  Now with f2.
f2 = ff(1)*P(1,:) + ff(2)*P(2,:);
plot(k,f,'k',k,f1,'r',k,f2,'g')
% which is pretty good.  How about w/ f3.
f3 = ff(1)*P(1,:) + ff(2)*P(2,:) + ff(3)*P(3,:);
plot(k,f,'k',k,f1,'r',k,f2,'g',k,f3,'b')
% and finally f4 (which is a 3rd order exp)
f4 = ff(1)*P(1,:) + ff(2)*P(2,:) + ff(3)*P(3,:) + ff(4)*P(4,:);;
plot(k,f,'k',k,f1,'r',k,f2,'g',k,f3,'b',k,f4,'c.')
legend('f','0','1','2','3',0)